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Abstract 



We use the Bogoliubov theory of atoms in an optical lattice to study 
the approach to the Mott-insulator transition. We derive an explicit ex- 
pression for the superfluid density based on the rigidity of the system 
under phase variations. This enables us to explore the connection be- 
tween the quantum depletion of the condensate and the quasi-momentum 
distribution on the one hand and the superfluid fraction on the other. The 
approach to the insulator phase may be characterized through the filling of 
i-^J ■ the band by quantum depletion, which should be directly observable via 

C"j ' the matter wave interference patterns. We complement these findings 

© , by self-consistent Hartree-Fock-Bogoliubov-Popov calculations for one- 

O ■ dimensional lattices including the effects of a parabolic trapping potential. 

> : 

X 1 1 Introduction 



Spectacular progress has been made in experimental studies of atoms loaded 
into an optical lattice in the region of the Mott superfluid insulator quantum 
phase transition ^[5]. In this article, we shall discuss the superfluid density 
and the quasi-momentum distribution, which is directly related to the matter- 
wave interference patterns that can be observed in such experiments. To do this 
we use the Bogoliubov method [3], as developed for use in optical lattices @]. 
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In a previous paper |5| we used this method to produce results for squeezing 
that are consistent with those of other approaches previously reported in the 
literature jS] \7\ |S1 E] • In this paper we shall show how it can be used to predict 
the decrease in the superfluid fraction and the corresponding variations in the 
matter wave interference fringes that should be directly observable in future 
experiments. This extends our previous studies based on exact calculation for 
small one-dimensional systems |10| into the experimentally relevant regime of 
lattice sizes and particle numbers. 

We first introduce the Bose-Hubbard Hamiltonian for atoms in an optical 
lattice QT] . We then describe briefly our method for determining the superfluid 
fraction based on the rigidity of the system under a twist of the condensate 
phase ^21- Using a perturbative formulation analogous to the Drude weight 
|13| the Bogoliubov approximation gives us a particularly direct way of finding 
this quantity. It also gives a simple picture of how superfluidity is suppressed by 
quantum depletion of the condensate. We shall compare the results for various 
quantities calculated using the Bogoliubov approximation with exact numerical 
calculations for the case of modest numbers of atoms and lattice sites |l(Jj . 

2 The Bose-Hubbard Model and Superfluidity 

The Bose-Hubbard Hamiltonian for atoms in an one-dimensional optical lattice 
with / sites has the form [TT] : 

11 V 1 

H = ^2ni€i - J ^2(a\ +1 ai + a\a i+1 ) + — - 1) . (1) 

i—1 i—1 i—1 

Here J represents the coupling between adjacent lattice sites due to tunneling 
and V is the strength of repulsion between atoms on the same site. The non- 
interacting energy of the atoms on each site, e*, will have some variation that is 
typically smooth on the scale of the condensate. We shall consider below both, 
the case where this is a constant, as well as the extension to the case where 
it varies in a trapped condensate. This Bose-Hubbard Hamiltonian should be 
an appropriate model when the loading process produces atoms in the lowest 
vibrational state of each well, with a chemical potential smaller than the distance 
to the first vibrationally excited state. This is known to be possible from the 
results of recent experiments [Tl l2l I14|. 

The concept of superfluidity is closely related to the existence of a condensate 
in the interacting many-body system. Formally, the one-body density matrix 
p^ 1 ) (x, x') has to have exactly one macroscopic eigenvalue which defines the 
number of particles in the condensate; the corresponding eigenvector describes 
the condensate wave function <fio (x) — e 1 ®^ \<j>Q (x)\. A spatially varying con- 
densate phase, (x) , is associated with a velocity field for the condensate by 



This irrotational velocity field is identified with the velocity of the superfluid 
flow, v s (x) = Vq (x) (dS],[Ill) and enables us to derive an expression for the 
superfluid fraction, f s . Consider a system with a finite linear dimension, L, 
in the e\ -direction and a ground-state energy, Eq, calculated with periodic 
boundary conditions. Now we impose a linear phase variation, (x) = 9xi/L 
with a total twist angle 9 over the length of the system in the ei-direction. 
The resulting ground-state energy, Eg will depend on the phase twist. For very 
small twist angles, 9 <C 7r, the energy difference, Eg — E , can be attributed to 
the kinetic energy, T s , of the superflow generated by the phase gradient. Thus, 

E e -E =T s = ^mNf s iP s , (3) 

where m is the mass of a single particle and N is the total number of particles 
so that mNf s is the total mass of the superfluid component. Replacing the 
superfluid velocity, v s with the phase gradient according to Eq. © leads to a 
fundamental relation for the superfluid fraction 



_ 2m L 2 Eg — Eq _ 1 Eg- Eq 

/s h 2 N 9 2 N J(A9) 2 ' 



where the second equality applies to a lattice system on which a linear phase 
variation has been imposed. Here the distance between sites is a, the phase 
variation over this distance is A9, and the number of sites is /. In this case, 
J=h 2 /(2ma 2 ). 

Technically the phase variation can be imposed through so-called twisted 
boundary conditions 12] . In the context of the discrete Bose- Hubbard model it 
is, however, more convenient to map the phase variation by means of a unitary 
transformation onto the Hamiltonian. The resulting "twisted" Hamiltonian 

11 V 1 

Hg = ^n l£l - J^(e- lAe ^ +1 a, + e^ 6 a\a l+1 ) + -J^n^i - 1) (5) 

i—l i— 1 i— 1 

exhibits additional phase factors e ±lAe — the so-called Peierls phase factors 
in the hopping term |17l 118) . These phase factors show that the twist is 
equivalent to the imposition of an acceleration on the lattice for a finite time. It 
is interesting to note that the present experiments enable us to make a specific 
connection between the formal and operational aspects of the system. 

We calculate the change in energy Eg — Eq under the assumption that the 
phase change A9 is small so that we can write: 

e~ iAe ~ 1 - iA9 ~ ^(A9) 2 . (6) 

Using this expansion the twisted Hamiltonian JSJ) takes the following form: 

He~H + A9J-^(A9) 2 f = H + H pcrt , (7) 
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where we retain terms up to second order in AO. The current operator J (N.B. 
that the physical current is given by this expression multiplied by t) and the 
hopping operator T are given by: 

J = iJ J2(aj +1 a 2 - a\a l+1 ) (8) 

i=l 
I 

f = -J^2(a\ +l di + ala i+1 ) . (9) 
i=i 

The change in the energy Eg — Eq due to the imposed phase twist can now be 
evaluated in second order perturbation theory 

E - E = AE {1) + A£ (2) . (10) 

The first order contribution to the energy change is proportional to the expec- 
tation value of the hopping operator 

AEW = (* |# port |*o) - ~^(A0) 2 (^ a \f |* ) . (11) 

Here |\&o) is the ground state of the original Bose-Hubbard Hamiltonian JJJ. 
The second order term is related to the matrix elements of the current operator 
involving the excited states \^f v ) (v = 1, 2, ...) of the original Hamiltonian 

±E (2)_ y l(^|gp CT t|* )| 2 _ (A6 n2 V l(^|j|*0)| 2 fl2 x 

^ E„-E { ' ^ Ev-Eq ' [ 1 

Thus we obtain for the energy change up to second order in AO 

|(*,|i|*o^ 12 



E e -E = (Afl) a (- 1(^1*0) -E e e )=im 2 D, 
D . (1 3, 

The quantity D, defined above, is formally equivalent to the Drude weight used 
to specify the DC conductivity of charged fermionic systems ^3] ■ The superffuid 
fraction is then given by the contribution of both the first and second order term: 

f - fM - f( 2 )- 
J s J s J s ' 

/ s (1) = -^j{(^\f\^o)), (14) 

m = 1 / v I(^I*q)I 2 \ 

/s NJ\2-i E v -E r 

Here N is the number of atoms in the lattice. In general both, the first and the 
second order term contribute. For a translationally invariant lattice the second 
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term vanishes (as is going to be shown latter) in the Bogoliubov limit that we 
shall use in this study. However, in exact calculations and in the Bogoliubov 
approximation for an inhomogeneous lattice the second order term plays a role. 

We can further understand this approach to the superfluid density by cal- 
culating the flow that is produced by the application of the phase twist. To 
do this we work out the expectation value of the current operator expressed in 
terms of the twisted variables: 

Jg = uY,{e- iA6 a\ +1 k ~ e iAe ala l+1 ) . (15) 
i=l 

We expand this to find the lowest order contributions, i.e.: 

i 

jg~j + JAe^ial^ + a\a l+1 ) = J — TAB . (16) 

i=l 

We use first order perturbation theory on the wave function to obtain the fol- 
lowing expression: 

(¥(A0)|J fl |¥(A0)) = 2A9( - ^(ttolTltto) - E ^ 

= 2NJf s A6. (18) 

If we note that the kinetic energy for a small quasi-momentum q on a lattice 
is given by Jq 2 a 2 , we can define the effective mass as m* = jj^t • Here the 
quasi- momenta are given by q = j^- j with j = 1, ...,(/ — 1) and lattice spacing 
a. Thus, the physical current, Eq. (|18f) multiplied by j-, can be expressed as: 

(*(A0)|J fl |*(A0)) = Nf s A6^ . (19) 

ma 

This is the total flux and we need to divide I to get the flux density, i.e. 



1 (*(a*m i*(a*)> - {*M\r N & 



I \m*a J \ al 



(20) 



So we see that the Drude formulation of the superfluid fraction (|14f> gives an 
intuitively satisfying expression for the amount of flowing superfluid. 

3 The Bogoliubov Approximation to the Bose- 
Hubbard Hamiltonian 

We use the Bogoliubov approximation for the Bose-Hubbard model in the limit 
that quantum fluctuations, or equivalently depletion of the condensate, is not 
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too great. In the limit that the quantum depletion can be completely ignored, 
we can replace the creation and annihilation operators, a\ and hi, on each site 
with a c- number, Zj. This leads to a set of coupled non-linear Schrodinger, i.e. 
Gross-Pitaevskii (GP), equations for these amplitudes |19j : 

ihd t Zi = -J(z i+ i + Zi-i) + VziZ*Zi . (21) 

This equation can be used to study the properties of the condensate loaded 
into the lattice when the tunneling kinetic energy is large enough compared 
to the interaction energy though small enough for the one-band Bose-Hubbard 
model to be valid. We then include the quantum fluctuations in our description 
of the system using the Bogoliubov approximation, where we suppose that we 
can write the full annihilation operator in terms of the c-number part and a 
fluctuation operator thus: 

a i = (z i + S i )e- i ^ . (22) 

This form will be useful when we are looking at the properties of a time- 
independent or adiabatic ground state. In using this method we are assuming 
that the fluctuation part is small. The Bogoliubov method gives us expres- 
sions for the averages of the squares of the fluctuation operator and allows us 
to determine whether this assumption is valid. We shall examine its validity 
by comparing the results for various physical quantities with exact numerical 
calculations based on the Bose-Hubbard Hamiltonian. 



3.1 Bogoliubov theory for the translationally invariant lat- 
tice 

The ground state solution for the translationally invariant lattice gives the eigen- 
value: 

H = n a V-2J, (23) 

where 

n = N/I (24) 

is the mean number of atoms on each site of the lattice. We take N to be the 
total number of atoms and / to be the number of sites in the one-dimensional 
lattice. 

The Bogoliubov equations for the lattice have the following form: 

md t 5i = (2n V - (1% - J(S i+l + k-x) + noVSj . (25) 

This is solved by constructing quasi-particlcs for the lattice which diagonalizc 
the Hamiltonian 4 , i.e 

Si = ^^[vP&qe**"-"**) - v q * a\ e-' l{qia -^] (26) 
Sj = ^ u9 * A l e-*"™-"^ - v q a q e^ qia -^ t) ] , (27) 
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where a is the lattice spacing. The quasi-particle operators obey the usual Bosc 
commutation relations: 

and have the following expectation values at some temperature T: 

{a\a ql ) = 5 qq ,[exp(huj q /k b T) - (29) 

We then find the following equations for the excitation amplitudes and fre- 
quencies, 



hu>„ 



-htu q v q 



^ + 4/ sin 2 (y) u q -n Vv q , 
n V + 4Jsin 2 (y) v q -n Vu q . 

Thus, the expressions for the u q and v q yield: 

, g,2 _ K (Q) + n V + huJ 1 



„?|2 



2hui q 
K(q) + n V — huj q 
2huj„ 



where the phonon excitation frequencies are given by: 

htv q = ^K(q)[2n V + K(q)} 
K(q) = 4 J sin 2 (^) . 



(30) 
(31) 

(32) 
(33) 



(34) 
(35) 



3.2 Expressions for the number superfiuid fraction in the 
translationally invariant lattice 

Having obtained the expressions for the excitations we can now determine the 
superfiuid fraction. The quantity we need to calculate is just the first order 
term of the Drude expression (|14p> . because the second order term vanishes in 
the Bogoliubov limit due to the translational invariance of the lattice [see Eq. 
JUL i-e. 



1 1 1 



(36) 



i=i 



In the Bogoliubov approximation this has the form 
i 

■is 



L = 7^I](*o|(^ + l+^+l)(^ + ^) + (^ + ^)(^+l+^+l)l*o) 
i=l 

i- £> |2z? + SlJi + $$ +1 |* ) • (37) 



2N 
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We can now express the fluctuation operators, Eqs. iffifl) and J22J), in terms of 
the quasi-particle operators that diagonalize the quadratic Hamiltonian. This 
leads to the following expression for the superfluid fraction at finite temperature: 



J 2N 



i=l 

X 



=1 q 

Qt/aJ, e-' iq ' la -v q 'a q , e iql 
q' 

(38) 



-w> _ „«a e iqia ] 



i—1 q 

x J^K'cv e lq ' {l+1)a - v q 'a\, e ~ iq ' {l+1)a ] 
<?' 

and we find in the zero temperature limit of a translationally invariant lattice: 

I 



fs v 



z 2 + 



l^|^| 2 cosM] . (39) 



Here the summation runs over all quasi- momenta q = ^Ljwithj = 1, ...,(/— 1) 
and we have called z the value of all Zi in a translationally invariant system. 
This shows that in the limit of zero lattice spacing (while keeping q finite) the 
superfluid fraction is unity as we have the normalization condition: 

Iz 2 + J2\v q \ 2 = N . (40) 

These expressions give a direct insight into the change of the superfluid fraction 
as atoms are pushed out of the condensate due to interactions. In Eq. l|3"9f 
the sum involving the Bogoliubov amplitudes v q characterizes the difference 
between the condensate fraction, which is given by the first term, and the su- 
perfluid fraction. For weak interactions and a small depletion, which fills only 
the lower quarter of the band where the cos(ga) term has a positive sign, the 
superfluid fraction is larger than the condensate fraction. Thus the depletion 
of the condensate has initially little effect on superfluidity When the depleted 
population spreads into the central part of the band, where the cos(qa) term 
has a negative sign, the superfluid fraction is reduced and might even become 
smaller than the condensate fraction. Finally, the population in the upper quar- 
ter of the band again produces a positive contribution to the superflow. In a 
sense the interactions are playing a role akin to Fermi exclusion "pressure" in 
the case of electron flow in a band. This, however can lead to perfect filling and 
cancellation of the flow. In the case of our Bogoliubov description we can only 
see reduction of the flow, not a perfect switching off of the superfluid. This hap- 
pens in the Mott insulator state, which cannot be described by the Bogoliubov 
approximation. 
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In the next section we outline the version of the Bogoliubov theory that 
should be best suited to treating these problems, i.e self-consistent Bogoliubov 
theory. 



4 Self-consistent HFB-Popov theory 

In this section we explore the limits of validity of the simplest zero tempera- 
ture self-consistent Bogoliubov theory, a simplified version of the Hartree-Fock- 
Bogoliubov approximation originally introduced by Popov [2U]. The HFB-Popov 
theory is an extension of the standard Bogoliubov approximation that takes into 
account the depletion of the condensate but neglects the anomalous average. 
As discussed in the previous section, taking into account the depletion of the 
condensate is important as the transition is approached because the depleted 
population causes the reduction of the superfluidity. Although the HFB-Popov 
aproach has the limitation that it doesn't take into account the full effect of the 
medium because it neglects the anomalous average [U, it can be considered a 
better theory for the elementary excitations than the full HFB due to the fact 
that the theory is gapless and doesn't violate Goldstone's theorem. 

A derivation of the Bogoliubov equations for the quasiparticle amplitudes in 
a lattice can be found for example in Ref. [Sj. These equations only take into 
account terms up to second order in the fluctuations. Including third and fourth 
order terms by treating them in a self-consistent mean field approximation ( |22|. 

) and neglecting anomalous average terms yields the following HFB-Popov 
equations: 

^ul + cHi = (2V(\z l \ 2 + fi l )~ f , + e t )u!~J(u! +1 +ul 1 )~VzM^) 
-hu q vl-c«z* = {2V{\z l \ 2 +n l )~^ + e l )vJ-J{v'l +1 +vU)-Vzf^2) 
ftZi = -J(z l+1 +z l _ 1 ) + (V(\z l \ 2 + 2fi,) + e t )z l , (43) 

Ek'l 2 ' ( 44 ) 



N = ^(N 2 +n 4 ), (45) 
i=i 

c« = V^^M-z^). (46) 

i 

Where ei is the energy offset at site i due to an external potential, {uf, vf} and 
ui q are respectively the quasiparticle amplitudes and energies, thus 



u\di q e-^ - vfa\e^ 1 , (47) 



hi is the depletion at site i, and N is the total number of particles. The param- 
eters c q ensure the {uf , if} solutions to the above equations with ui q ^ to be 
orthogonal to the condensate (Ref. \'2'6\). 
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By calculating the quasiparticle amplitudes and the condensate density it is 
possible to get information about most of the physical properties of the system. 
For example the superfluid fraction and the on site number fluctuations are 
given by: 

f s = / s (1) -/i 2) , (48) 




f (2) = ±_ \ i £m - u l v i+i)\ 2 + s / lEiOf+i^ - u i v i+i)\ 

TlOJq + TlUJq' qq 2flUJ q 

An? = W 2 X;K-«ri 2 - (49) 
1 

From the complete expression of the superfluid fraction, it can be seen explicitly 
how due to the translational invariance, the second order term vanishes in the 
homogeneous system. 



4.1 Translationally Invariant lattice 

For the translational invariant lattice we use the quasiparticle transformation 
given by Eqs. (|26|l and 1)27(1 . Under this transformation the self consistent 
equations can be written, generalizing the previous version, as: 



A 4 = \\z\' + 1 ^\vrjV-2J, (50) 
K(q) + \z\ 2 V + hu q 



2flU)q 



(51) 



a _ K(g) + \z\ 2 V-hu, q 



2n» q (52) 

N = I\z\ 2 + J2\v q \ 2 - (53) 
i 

Here the phonon excitation spectrum is given by: 

JUj q = ^/K(q)[2\z\ 2 V + K(q)}, (54) 

and K(q) is given by Eq. (|35[1 ■ Again we omit the subscript in the amplitudes 
Zi because they have the same value at all lattice sites. Notice also that due to 
translational invariance, the c q coefficients vanish. 

In the homogeneous system, the form of the HFB-Popov equations for the 
quasiparticle amplitudes and energies is very close to the standard Bogoliubov 
form. We do, however, have to replace no = N/I by the condensate amplitude 
| z | 2 which must take into account the depletion of the condensate. We solve 
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for the condensate amplitude as a function of the external parameters J, V, N 
and I by inserting Eq. I|52l) in Eq. (|53[l . Once \z\ 2 is determined, we use it to 
calculate the other expressions. 

In Fig. n we compare the number fluctuations on a lattice site, the con- 
densate fraction and the total and second order superfluid fraction determined 
from the exact solution of the Bose-Hubbard Hamiltonian to the self consistent 
HFB-Popov predictions as a function of the ratio V e // = V/J. The systems 
used for the comparisons have three wells, 1 = 3, and commensurate filling 
factors uq = 5, 10, 20 and 50. Similar results for the incommensurate case with 
N = 16,31,61, 151 are shown in Fig. |3 We were restricted to consider only 
three wells due to computational limitations. The size of the matrix needed in 
the exact solution for N atoms and / wells scales as jjTTfErjr ■ However, if the 
HFB-Popov approach works well for these small systems we expect it to provide 
a good description of the larger systems prepared in the lab. 

Because the second order term of the superfluid fraction (second term of 
Ea ll4f) vanishes in the HFB-Popov approach (see Eq. 148(1 . we only expect the 
self consistent HFB-Popov theory to give a good description of the superfluid 
fraction in the region where the second order term is extremely small, provided 
it predicts accurately the first order term. This is exactly what is observed in the 
plots. When the second order term starts to grow, typically above 0.5V£fj, the 
HFB-Popov theory starts to fail. An estimate of is shown by a vertical line 
in some of the figures. This was obtained by using the second order perturbative 
approach presented in Ref. 0]. With increasing filling factor the critical value 
is shifted towards larger values of the interaction strength, and the region in 
which the HFB-Popov theory is accurate gets larger. It is interesting to note 
that the number fluctuations predicted by the theory are accurate in a greater 
range than the other physical quantities shown. Its predictions of squeezing 
agree very well with the exact solutions right up to the point where the number 
fluctuations become less than unity. 

For the cases with non-commensurate fillings depicted in Fig. [2J the agree- 
ment is significantly better for all quantities. This is not surprising because 
when the filling is not commensurate there is always a superfluid present and 
the Mott transition doesn't occur. As can be seen in the plots for these cases 
the second order term is always very small. 

4.2 Inhomogeneous lattice 

In this section we consider the experimentally relevant case when there is an 
external magnetic confinement in addition to the lattice potential. In this situ- 
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Figure 1 : Comparisons of the exact solution (dotted line) and HFB-Popov (solid 
line) as a function of V e // = V/J, for a system with 1 = 3 and filling factors 
n D = 5, 10, 20, and 50. Top: number fluctuations, middle: condensate fraction, 
bottom: superfluid fraction / s . The exact second order term (dashed line) of 

the superfluid fraction, / s is also shown in these plots. The vertical line shown 
in some plots is an estimation of V^f ■ 



ation, the self consistent HFB-Popov equations take the form: 



flUJ q U q i +C q Zi 



-TVjJqVl 



hi 
N 



= (2V(\z l \ 2 +n i ) - M + Sli 2 )ul- J(u q i+1 +uU) - Vz 2 v%55) 
= (2V(\z i \ 2 +h i ) - M + Qi 2 )v g - J(v q i+1 +vU) - Vz?u$6) 



-J(z i+1 +z i -. 1 )+ (V(\z t \ 2 



2h, 



i 



* q 



ZiV q ) 



(57) 
(58) 

(59) 

(60) 



where f2 = ^muj 2 a 2 , with m the mass of the atoms, u> the trap frequency, and a 
the lattice spacing. The site indices i are chosen such that i = corresponds to 
the center of the trap. Again the c q are introduced to ensure the orthogonality of 
the excitations to the condensate • We have solved the HFB-Popov equations 
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Figure 2: Comparisons of the exact solution (dotted line) and HFB-Popov 
(solid line) for a system with 1 — 3 and non commensurate filling factors 
TV = 16,31,61,151 as a function of V e ff — V/J. Top: number fluctuations, 
middle: condensate fraction, bottom: superfluid fraction / s . In these plots the 
exact second order term of the superfluid fraction is also shown with a dashed 
line. 



for this system by an iterative procedure, similar to the one followed in Ref. 
|24|. Each cycle of the iteration consists of two steps. In the first step we 
solve Eq. fBTjl subject to the constraint Eq. (fB^Qf by using the hi obtained in 
the previous cycle. This generates new values for the z,-. In the second step 
we solve for {uf,vf} in Eqs. ltB^|l using the hi from the previous cycle and the 
newly generated Zi. The {uf, vf} are used then to update hi. Because the HFB- 
Popov is gapless, it is possible to keep the orthogonality of the excitations to 
the condensate by solving Eqs. (|5^|l with the c q set to zero but removing in each 
cycle the projection of the calculated {uf,Vi} amplitudes onto the condensate. 
Convergence is reached when the change in \ni\ 2 from one cycle to the next 
is smaller than a specified tolerance. 

The parameters chosen for the numerical calculations were f2 = 0.0015-Er, 
with En the one photon recoil energy, which for the case of a Rubidium con- 
densate corresponds to a trap frequency of approximately 90 Hz. We used a 
total number of 1000 atoms, N = 1000, and set VN = 1.0-Er. J was varied to 
achieve a range of V e ff = V/J between 0.01 and 312. The range was chosen 
based on a local mean field approach 0], which for our parameters estimates 
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the transition region between V e ff » 640 (at the center where the local filling 
factor is approximately 80) and V e ff « 12 (at the wings). 

The results of the numerical calculations are summarized in Figs. |21to|SJ In 
Fig. |3 we plot the evolution of the density profile (black boxes), the condensate 
population (triangles) and the on-site depletion (empty diamonds) as V e ff is 
increased. In the plots we also show, for comparison purposes, the ground state 
density profile for J — (empty boxes). This has the advantage that can be 
calculated exactly from the Hamiltonian. In general we observe the reduction 
of the condensate population and thus the increment of the depletion as the 
interaction strength is increased. When the system is in the superfluid regime 
most of the atoms are in the condensate but as J is decreased the depletion of 
the condensate becomes very important. 

For the chosen parameters, the density profile has a parabolic shape reflecting 
the confining potential. By comparing the evolution of the density as J is 
decreased with the exact solution at J = 0, we can crudely estimate the validity 
of the HFB-Popov calculations. The density evolves from a Gaussian type (see 
plots for V e ff = 0.01 and 0.09) with smooth edges towards a Thomas-Fermi 
profile with sharp edges adjusting its shape to the J = profile. We can 
appreciate that around V e ff = 3 both profiles are almost equal. For lower 
values of J the HFB-Popov density starts to differ from the J = one, even 
though the system is closer to the J = limit. We can say that beyond this 
point higher order correlations, neglected by the theory, begin to be important. 
The departure of the HFB-Popov density profile from the J = one as J 
is decreased begins at the edges (see panel corresponding to V e ff — 11 and 
100). This is something expected if we look at the on-site depletion. For such 
values of Ve.fi the local depletion at the wings corresponds to a considerable 
percentage of the condensate populations, and thus the validity of the HFB- 
Popov assumptions starts to be dubious. The homogeneous results shown in the 
previous section corroborate our present statements for the confined system. For 
the smallest filling factor (see FigQ the differences between the homogeneous 
HFB-Popov calculations and the exact solutions become important for values 
of V e ff greater than 20. For higher values of V e ff, see plot for V e ff = 312, 
the HFB-Popov density predictions differs from the J = solution even at the 
central wells. At this point the failure of the method is clear and a fully quantal 
method is required. 

The HFB-Popov quasiparticle spectrum is shown in Fig0] It can be observed 
how the lower energy eigenvalues evolve from a linear non degenerated spectrum 
to an almost degenerated one as J is decreased. It is worth it to mention that 
the small energy difference between the ground and first excited states for high 
values of V e f / makes the numerical solution very unstable in the sense that it 
was very easy to jump to an excited state when solving for the condensate wave 
function. The decrement in the energy spacing predicted by the HFB-Popov 
theory as the system approaches the transition is very useful to keep in mind 
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for the experimental realization of the Mott transition. As the optical lattice 
depth is ramped up the adiabaticity criteria is harder to fulfill. 

In Fig. [S] we plot the results for the number fluctuations found numeri- 
cally using the inhomogeneous HFB-Popov approach. The number fluctuations 
profile reflects the condensate profile. We also show the number fluctuations 
evaluated by using a local density approximation (empty boxes). The latter 
was calculated by substituting in the number fluctuations expression (Eq. I49fl 
the {u q ,v q } amplitudes found for the homogeneous system (Eqs. ISH and . 
but replacing the condensate density in each lattice site by the one found nu- 
merically for the trapped system (see Fig. The complete agreement between 
the two approaches justifies the validity of the local density approximation for 
the estimations of local quantities in confined systems. Based on this agreement 
and the results for the homogeneous system shown in the previous section, we 
expect that the inhomogeneous HFB-Popov results for squeezing also agrees 
with the exact solution right up to the transition. 



4.3 Quasi-momentum distribution in the inhomogeneous 



The quasi-momentum distribution of the atoms released from the lattice is im- 
portant because it is one of the most easily accessible quantities to the experi- 
ments. The quasi-momentum distribution function n q is defined as |10| 



where the quasi-momentum q can assume discrete values which are integer mul- 
tiples of j^, a is the lattice spacing. In Fig. [B]we present the quasi-momentum 
distribution for the same parameters used in the previous section. The distribu- 
tion for the two lowest values of V e f / corresponds to the one that characterizes 
an uncorrelated supcrfluid phase with a narrow peak at small quasi-momcnta. 
As the hopping rate is decreased we observe that the sharpness of the central 
peak decreases and the distribution extends towards large quasi-momenta. It is 
interesting to note the appearance of a small peak between q = 0.5 and 1 which 
is most noticeable for the V e ff = 3 case. This agrees with the results found 
in |25j where they solve numerically the Bose-Hubbard Hamiltonian by using 
Monte Carlo simulations. We attribute the origin of the small peak to the deple- 
tion of the condensate at the wings. For the parameters when the small peak is 
present, the most important contribution to the quasi-momentum distribution 
still comes from the condensate atoms. The step function like shape of the con- 
densate profile causes an oscillatory |sin(a;)/a;| shape of the quasi-momentum 
distribution. As the lattice depth is increased the hopping becomes energeti- 
cally costly, the long-range order starts to decrease and the Fourier spectrum 
becomes broader. 



lattice 




(61) 
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In Fig. we plot the first order on site superfluid fraction f\ ' which was 
defined in Ea. (|48|l . The curves corresponding to V e ff = 0.01 — 11, which are 
in the regime where the HFB-Popov is expected to be valid, depict how as 
V e f f is increased the superfluid profile decreases faster at the wings and at the 
center but no major change is observed in the middle section. The evolution 
of the on-site superfluidity as the interaction strength is increased, exhibiting 
a domain localized decrement instead of a global one, is in agreement with the 
development of uncompressible regions surrounded by superfluid rings predicted 
for trapped systems |26| as the transition is approached. 

Fig. |H| shows the first order and total superfluid fraction and also the second 
order superfluid fraction as a function of V e ff. Different from the translation- 
ally invariant case, the second order contribution calculated in the HFB-Popov 
approach doesn't vanish for the inhomogeneous system. The rapid decrement of 
the superfluid fraction observed after Log(V e ff) ~ 1.2 is a signature that the 
system is entering a highly correlated regime. Beyond this point higher order 
correlations neglected in the HFB-Popov approach become crucial and a more 
accurate approach is required. 

The Mott transition is a quantum phase transition and as all critical phe- 
nomena its behavior depends strongly on the dimensionality of the system. In 
the present analysis, due to computational limitations, we considered one di- 
mensional systems. Experimentally, the Mott transition has been achieved 2 
in a 3 dimensional lattice with filling factors between 1 and 3. Even though the 
HFB-Popov approach fails to describe the strong coupling regime for the one 
dimensional systems we considered in the present paper, we showed how the 
method is incredibly powerful in describing most of its characteristic features as 
they are driven from the superfluid regime towards the transition. We expect 
the HFB-Popov method to give a better description of the transition as the 
dimensionality of the system is increased and therefore to be a good model in 
an experimental situation. 

As shown in previous studies , |2H] the Mott transition in a d-dimensional 
homogeneous system has two different critical behaviors: one (d+1) XY- like, 
for systems with fixed integer density as the interaction strength is changed, 
and one mean field-like exhibited when the transition is induced by changing 
the density. Different from the homogeneous case where the Mott transition 
is characterized by the global offset of the superfluidity, for confined systems, 
commensuration is only well defined locally. The inhomogeneity introduced by 
the confined potential allows the existence of extended Mott domains (above a 
critical interaction strength) surrounded by superfluid ones [26], thus the total 
superfluid fraction doesn't vanishes in the Mott regime. This issue, together 
with the fact that the finite length scale introduced by the trap suppresses the 
long wave fluctuations which are responsible for destroying the mean field [23I 1 . 
make us believe the critical behavior in confined systems to be more mean-field 
like. Because the critical dimension for the latter type of transition is two |27|. 

x One obvious consequence of this is that BEC is possible in one and two dimensions in a 
trap whereas in the homogeneous, thermodynamic limit it can not occur in fewer than three 
dimensions 
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|2*8). we expect that for trapped systems in d = 3, the range of validity of the 
HFB-Popov extends closer to the transition. 

5 Summary 

We have developed in this article a Bogoliubov method for describing the ap- 
proach of a condensate loaded in an optical lattice towards the Mott transition. 
We have shown that this method can be used to predict the relevant physical 
quantities over a useful range. We have also shown how it gives a powerful 
insight into the way quantum depletion reduces the long range order and the 
superfiuid fraction. 
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Figure 3: Condensate density (triangles), total density (filled boxes) and local 
depletion (empty diamonds) as a function of the lattice site for different values of 
V e ff. Although these quantities are defined only at the discrete lattice sites we 
join them to help visualization. The empty boxes represent the exact solution 
for the case J=0. 
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Figure 4: Quasiparticle spectrum in ascending order predicted by the HFB- 
Popov theory for different values of V e ff. Empty diamonds (V e /f — 0.01), stars 
(V e ff = 0.09), crosses (V e ff = 3), filled diamonds (V e ff = 11), empty boxes 
(V ef f = 100) and polygons (V ef f = 312) 
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Figure 5: Number fluctuations in the self consistent HFB-Popov approach as a 
function of lattice site for V e ff = 0.01 (boxes), V e ff = 0.09 (crosses), V e // = 3 
(circles), 1 ^// = H (triangles), V e // = 100 (stars) and V e ff = 312 (diamonds). 
The maximum value reached by the profile decreases as V e ff is increased. The 
empty boxes shown for each of the curves correspond to the number fluctua- 
tions predicted by the homogeneous HFB-Popov model using a local density 
approximation. 
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Figure 7: First order on-site superfluid fraction as a function of the lattice site for 
different values of V e ff. Filled boxes: V e ff. =0.01, empty boxes: 14//. = 0.09, 
empty diamonds: 14//. = 3, stars: 14//. = 11, crosses: 14//. = 100 and 
triangles: 14//. = 312. 
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Figure 8: Top panel: First order (boxes) and total (stars) superfluid fraction 
as a function of V e ff. ■ Bottom panel : second order superfluid fraction as a 
function of V e f /. All these quantities are calculated in the self consistent HFB- 
Popov approach. 
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Abstract 

We use the Bogoliubov theory of atoms in an optical lattice to study 
the approach to the Mott-insulator transition. We derive an explicit ex- 
pression for the superfluid density based on the rigidity of the system 
under phase variations. This enables us to explore the connection be- 
tween the quantum depletion of the condensate and the quasi-momentum 
distribution on the one hand and the superfluid fraction on the other. The 
approach to the insulator phase may be characterized through the filling of 
the band by quantum depletion, which should be directly observable via 
the matter wave interference patterns. We complement these findings 
by self-consistent Hartree-Fock-Bogoliubov-Popov calculations for one- 
dimensional lattices including the effects of a parabolic trapping potential. 



1 Introduction 

Spectacular progress has been made in experimental studies of atoms loaded 
into an optical lattice in the region of the Mott superfluid insulator quantum 
phase transition [1, 2]. In this article, we shall discuss the superfluid density 
and the quasi-momentum distribution, which is directly related to the matter- 
wave interference patterns that can be observed in such experiments. To do this 
we use the Bogoliubov method [3], as developed for use in optical lattices [4]. 
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In a previous paper [5] we used this method to produce results for squeezing 
that are consistent with those of other approaches previously reported in the 
literature [6, 7, 8, 9]. In this paper we shall show how it can be used to predict 
the decrease in the superfluid fraction and the corresponding variations in the 
matter wave interference fringes that should be directly observable in future 
experiments. This extends our previous studies based on exact calculation for 
small one-dimensional systems [10] into the experimentally relevant regime of 
lattice sizes and particle numbers. 

We first introduce the Bosc-Hubbard Hamiltonian for atoms in an optical 
lattice [11]. We then describe briefly our method for determining the superfluid 
fraction based on the rigidity of the system under a twist of the condensate 
phase [12]. Using a perturbative formulation analogous to the Drudc weight 
[13] the Bogoliubov approximation gives us a particularly direct way of finding 
this quantity. It also gives a simple picture of how superfluidity is suppressed by 
quantum depletion of the condensate. We shall compare the results for various 
quantities calculated using the Bogoliubov approximation with exact numerical 
calculations for the case of modest numbers of atoms and lattice sites [10]. 

2 The Bose-Hubbard Model and Superfluidity 

The Bosc-Hubbard Hamiltonian for atoms in an one-dimensional optical lattice 
with / sites has the form [11]: 

11 V 1 

H = ~ J ^2(4+1^ + ^+1) + — ^h t {h t - 1) . (1) 

i—1 i—1 i—1 

Here J represents the coupling between adjacent lattice sites due to tunneling 
and V is the strength of repulsion between atoms on the same site. The non- 
interacting energy of the atoms on each site, ej, will have some variation that is 
typically smooth on the scale of the condensate. We shall consider below both, 
the case where this is a constant, as well as the extension to the case where 
it varies in a trapped condensate. This Bose-Hubbard Hamiltonian should be 
an appropriate model when the loading process produces atoms in the lowest 
vibrational state of each well, with a chemical potential smaller than the distance 
to the first vibrationally excited state. This is known to be possible from the 
results of recent experiments [1, 2, 14]. 

The concept of superfluidity is closely related to the existence of a condensate 
in the interacting many-body system. Formally, the one-body density matrix 
(x, x') has to have exactly one macroscopic eigenvalue which defines the 
number of particles in the condensate; the corresponding eigenvector describes 
the condensate wave function O (x) = e 1 ®^ \cj>o (x)\. A spatially varying con- 
densate phase, 6 (x), is associated with a velocity field for the condensate by 



This irrotational velocity field is identified with the velocity of the superfluid 
flow, v s (x) = v (x) ([15], [16]) and enables us to derive an expression for the 
superfluid fraction, f s . Consider a system with a finite linear dimension, L, 
in the e*i -direction and a ground-state energy, Eq, calculated with periodic 
boundary conditions. Now we impose a linear phase variation, O (x) = 9xi/L 
with a total twist angle 9 over the length of the system in the e\-direction. 
The resulting ground-state energy, Eg will depend on the phase twist. For very 
small twist angles, 9 <C n, the energy difference, Eg — E , can be attributed to 
the kinetic energy, T s , of the superflow generated by the phase gradient. Thus, 

E e -E =T s = ^mNf a v 2 s , (3) 

where m is the mass of a single particle and N is the total number of particles 
so that mNf s is the total mass of the superfluid component. Replacing the 
superfluid velocity, v s with the phase gradient according to Eq. (2) leads to a 
fundamental relation for the superfluid fraction 



_2mL 2 Eg-E _ 1 Eg - E 
h h 2 N 9 2 Nj^ef 



where the second equality applies to a lattice system on which a linear phase 
variation has been imposed. Here the distance between sites is a, the phase 
variation over this distance is A9, and the number of sites is /. In this case, 
J=h 2 /(2ma 2 ). 

Technically the phase variation can be imposed through so-called twisted 
boundary conditions [12]. In the context of the discrete Bose-Hubbard model it 
is, however, more convenient to map the phase variation by means of a unitary 
transformation onto the Hamiltonian. The resulting "twisted" Hamiltonian 

11 V 1 

He = - J 5>~ iA< ^+i^ + eiAe «I«*+i) + -^-$>i(rii - 1) (5) 

i— 1 i— 1 i— 1 

exhibits additional phase factors e ±lAe — the so-called Peierls phase factors 
in the hopping term [17, 18]. These phase factors show that the twist is 
equivalent to the imposition of an acceleration on the lattice for a finite time. It 
is interesting to note that the present experiments enable us to make a specific 
connection between the formal and operational aspects of the system. 

We calculate the change in energy Eg — Eq under the assumption that the 
phase change A9 is small so that we can write: 

e - iAe ~l-iA#-i(A0) 2 . (6) 

Using this expansion the twisted Hamiltonian (5) takes the following form: 

He ~H + A9J - ^(A9) 2 f = H + H pert , (7) 
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where we retain terms up to second order in A9. The current operator J (N.B. 
that the physical current is given by this expression multiplied by j-) and the 
hopping operator T are given by: 

J = U^ j {a\ +1 a i -a\a i+1 ) (8) 

i=l 

f = -jJ2(4+i^ + ala l+1 ) . (9) 

i=l 

The change in the energy Eg — Eq due to the imposed phase twist can now be 
evaluated in second order perturbation theory 

E -E o = A£ (1) + A£ (2) . (10) 

The first order contribution to the energy change is proportional to the expec- 
tation value of the hopping operator 

A£« = (*„ Import |* > = -i(A0) 2 (vI/ o |f |* ) . (11) 

Here |^o) is the ground state of the original Bose-Hubbard Hamiltonian (1). 
The second order term is related to the matrix elements of the current operator 
involving the excited states 1^^) (y — 1, 2, ...) of the original Hamiltonian 

AS(2) = _ y KWc^o)! 2 = _ {Mf y VMJ^l 

^-f E v - Eq *-f E v - E 

Thus we obtain for the energy change up to second order in A0 

|(*,|J|*o x|2 



Eg — Eq = (A9f ( - i<g„|f |*„) - Y: ) = IWfD, 

D . ^.l {M -^<^). (13) 

The quantity D, defined above, is formally equivalent to the Drude weight used 
to specify the DC conductivity of charged fermionic systems [13]. The superfiuid 
fraction is then given by the contribution of both the first and second order term: 

f - /"(I) _ f(2). 

/ s (1) = ( 14 ) 

Js NJ\^f n E v -E )' 

Here N is the number of atoms in the lattice. In general both, the first and the 
second order term contribute. For a translationally invariant lattice the second 
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term vanishes (as is going to be shown latter) in the Bogoliubov limit that we 
shall use in this study. However, in exact calculations and in the Bogoliubov 
approximation for an inhomogeneous lattice the second order term plays a role. 

We can further understand this approach to the superfluid density by cal- 
culating the flow that is produced by the application of the phase twist. To 
do this we work out the expectation value of the current operator expressed in 
terms of the twisted variables: 

i 

Jo = ijJ2(e- [A9 al+i^ ~ e iAe a\a l+1 ) . (15) 
We expand this to find the lowest order contributions, i.e.: 

j ~j + JAeJ2(a\ +1 ai + ojoi+i) = J — TA6 . (16) 

i=l 

We use first order perturbation theory on the wave function to obtain the fol- 
lowing expression: 

<*(A0)|J fl |*(A0)) = 2A6( l^olf^o) E K ^'^ )|2 ) 

= 2NJf s A6. (18) 

If we note that the kinetic energy for a small quasi-momentum q on a lattice 
is given by Jq 2 a 2 , we can define the effective mass as m* = jj^-- Here the 
quasi- momenta are given by q = -p- j with j = 1, ...,(/ — 1) and lattice spacing 
a. Thus, the physical current, Eq. (18) multiplied by j-, can be expressed as: 

(tf (A0)| J fl |*(A0)> = Nf s A8-^ . (19) 
This is the total flux and we need to divide I to get the flux density, i.e. 

= v s n s . (20) 

So we see that the Drude formulation of the superfluid fraction (14) gives an 
intuitively satisfying expression for the amount of flowing superfluid. 

3 The Bogoliubov Approximation to the Bose- 
Hubbard Hamiltonian 

We use the Bogoliubov approximation for the Bose-Hubbard model in the limit 
that quantum fluctuations, or equivalently depletion of the condensate, is not 
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too great. In the limit that the quantum depletion can be completely ignored, 
we can replace the creation and annihilation operators, a\ and fij, on each site 
with a c- number, z,. This leads to a set of coupled non-linear Schrodinger, i.e. 
Gross-Pitaevskii (GP), equations for these amplitudes [19]: 

ihd t z l = -J(z i+ i + Zi-i) + VziZ*Zi . (21) 

This equation can be used to study the properties of the condensate loaded 
into the lattice when the tunneling kinetic energy is large enough compared 
to the interaction energy though small enough for the one-band Bose-Hubbard 
model to be valid. We then include the quantum fluctuations in our description 
of the system using the Bogoliubov approximation, where we suppose that we 
can write the full annihilation operator in terms of the c-number part and a 
fluctuation operator thus: 

a* - {zi + 5i)e-^ . (22) 

This form will be useful when we are looking at the properties of a time- 
independent or adiabatic ground state. In using this method we are assuming 
that the fluctuation part is small. The Bogoliubov method gives us expres- 
sions for the averages of the squares of the fluctuation operator and allows us 
to determine whether this assumption is valid. We shall examine its validity 
by comparing the results for various physical quantities with exact numerical 
calculations based on the Bose-Hubbard Hamiltonian. 



3.1 Bogoliubov theory for the translationally invariant lat- 
tice 

The ground state solution for the translationally invariant lattice gives the eigen- 
value: 

H = n V - 2 J, (23) 

where 

n = N/I (24) 

is the mean number of atoms on each site of the lattice. We take N to be the 
total number of atoms and I to be the number of sites in the one-dimensional 
lattice. 

The Bogoliubov equations for the lattice have the following form: 

ihdA = (2n V - n)6i - J(5 i+1 + 8^) + n V5\ . (25) 

This is solved by constructing quasi-particles for the lattice which diagonalize 
the Hamiltonian [4], i.e 

& = -L VVdg e 1 ^-"**) - v q *a\ e -i(fl<a-««*)] (26) 
v q 

5} = 4= Y[u q *&l e -i(9*°-<".*) _ v ia q e i (9 fa -'"«*)] , (27) 
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where a is the lattice spacing. The quasi-particle operators obey the usual Bose 
commutation relations: 

[a q ,a\]=5 qq , (28) 
and have the following expectation values at some temperature T: 



(ajay) = <W [exp(ftw 9 /fc b T) - 1] 1 . 



(29) 



We then find the following equations for the excitation amplitudes and fre- 
quencies, 



hoj„ 



— hw q V q = 



n t/ + 4Jsin 2 (^) 
n V + 4Jsin 2 (f) 

Thus, the expressions for the u q and v q yield: 

i qi 2 K (q) + n oV + fe^g 
|U 1 - 2fuv q 

K(q) + n V - HuJq 

2htd„ 



u q - n Vv q , 
v q - n Vu q . 



„9|2 = 



where the phonon excitation frequencies are given by: 

hcoq = ^K{q)[2n V + K{q)} 
K(q) = 4Jsin 2 (f). 



(30) 
(31) 

(32) 
(33) 



(34) 
(35) 



3.2 Expressions for the number superfiuid fraction in the 
translationally invariant lattice 

Having obtained the expressions for the excitations we can now determine the 
superfiuid fraction. The quantity we need to calculate is just the first order 
term of the Drude expression (14), because the second order term vanishes in 
the Bogoliubov limit due to the translational invariance of the lattice [see Eq. 
(48)], i.e. 



1 1 1 



2N 



(36) 



i=l 



In the Bogoliubov approximation this has the form: 
1 1 

f> = ^S<*ol(*i+l + ^+l)(*i + «*) + (*i +«*)(*i+l + «i+l)l*o) 



i=l 



= ^E(*o|2z 2 + ^ +1 <5, + ^ +1 |*o) 



(37) 



i=l 
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We can now express the fluctuation operators, Eqs. (26) and (27), in terms of 
the quasi-particle operators that diagonalize the quadratic Hamiltonian. This 
leads to the following expression for the superfluid fraction at finite temperature: 



x > \w a', e iqia - v q a a , e lq m 



f> = i[E^+i(E[ u S ei,(,+1)a -^^ (,+1)a ] 

L 2—1 q 

5>*'4 e ~ iq 

+ jE(E[ u ^ e- iqia -v q a q e iqia ] 

i=l q 

x ^[u q 'a q , e lq ' {l+1)a - v q 'a\, e -V(.i+V°]} , (38) 

q' 

and we find in the zero temperature limit of a translationally invariant lattice: 

I 



Is N 



z 2 + 



±5><fcos(«a)] . (39) 
i 



Here the summation runs over all quasi- momenta q= j^j with j = 1 , . . . , (7 — 1) 
and we have called z the value of all z% in a translationally invariant system. 
This shows that in the limit of zero lattice spacing (while keeping q finite) the 
superfluid fraction is unity as we have the normalization condition: 

Iz 2 + ^\v q \ 2 = N . (40) 



These expressions give a direct insight into the change of the superfluid fraction 
as atoms are pushed out of the condensate due to interactions. In Eq. (39) 
the sum involving the Bogoliubov amplitudes v q characterizes the difference 
between the condensate fraction, which is given by the first term, and the su- 
perfluid fraction. For weak interactions and a small depletion, which fills only 
the lower quarter of the band where the cos(ga) term has a positive sign, the 
superfluid fraction is larger than the condensate fraction. Thus the depletion 
of the condensate has initially little effect on superfluidity. When the depleted 
population spreads into the central part of the band, where the cos(qa) term 
has a negative sign, the superfluid fraction is reduced and might even become 
smaller than the condensate fraction. Finally, the population in the upper quar- 
ter of the band again produces a positive contribution to the superflow. In a 
sense the interactions are playing a role akin to Fermi exclusion "pressure" in 
the case of electron flow in a band. This, however can lead to perfect filling and 
cancellation of the flow. In the case of our Bogoliubov description we can only 
see reduction of the flow, not a perfect switching off of the superfluid. This hap- 
pens in the Mott insulator state, which cannot be described by the Bogoliubov 
approximation. 
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In the next section we outline the version of the Bogoliubov theory that 
should be best suited to treating these problems, i.e self-consistent Bogoliubov 
theory. 



4 Self-consistent HFB-Popov theory 

In this section we explore the limits of validity of the simplest zero tempera- 
ture self-consistent Bogoliubov theory, a simplified version of the Hartree-Fock- 
Bogoliubov approximation originally introduced by Popov [20] . The HFB-Popov 
theory is an extension of the standard Bogoliubov approximation that takes into 
account the depletion of the condensate but neglects the anomalous average. 
As discussed in the previous section, taking into account the depletion of the 
condensate is important as the transition is approached because the depleted 
population causes the reduction of the superfluidity Although the HFB-Popov 
aproach has the limitation that it doesn't take into account the full effect of the 
medium because it neglects the anomalous average [21], it can be considered a 
better theory for the elementary excitations than the full HFB due to the fact 
that the theory is gapless and doesn't violate Goldstone's theorem. 

A derivation of the Bogoliubov equations for the quasiparticle amplitudes in 
a lattice can be found for example in Ref. [5]. These equations only take into 
account terms up to second order in the fluctuations. Including third and fourth 
order terms by treating them in a self-consistent mean field approximation ( [22] , 
[23] ) and neglecting anomalous average terms yields the following HFB-Popov 
equations: 

nu q u g i +c"z i = (2V(\z l \ 2 + f ll )- f i + e t )u^J(u! +1 +u!_ 1 )-Vzfvm) 
-nu q v1-ciz* = (2V(\z l \ 2 +h t )-^ + e l )v q - J{v1 +l +vU)- Vz?d$2) 
HZi = -J(z l+1 +z l - 1 ) + (V(\z t \ 2 + 2fi l ) + e t )z t , (43) 

EK 9 ! 2 . ( 44 ) 



n, 



N = ]T(N 2 +n 4 ), (45) 
= T/]>>| 2 (z>?-z^). (46) 

i 

Where is the energy offset at site i due to an external potential, {itf , vf} and 
uj q are respectively the quasiparticle amplitudes and energies, thus 

« < = 2«?a,e- i "« t - V ?*aJe iw « t , (47) 
<? 

hi is the depletion at site i, and N is the total number of particles. The param- 
eters c q ensure the {uf , v\ } solutions to the above equations with oj q ^ to be 
orthogonal to the condensate (Ref. [23] ) . 
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By calculating the quasiparticle amplitudes and the condensate density it is 
possible to get information about most of the physical properties of the system. 
For example the superfluid fraction and the on site number fluctuations are 
given by: 

fs = fP-fP, (48) 

/ s (i) = e^ 1)= ^eM+^)+e(«i+w]' 



i=i 



J l^lEii^vf - u y i+1 )\ 2 t x |£>? +1 <-^ +1 )' 2 



f(2) t y L^j \^+i u i Z ^iiwj I x 

N \jj hw q + huj ql qq ' 2huo q 



An? = N 2 £l«?-t\ ? T- (49) 
<i 

From the complete expression of the superfluid fraction, it can be seen explicitly 
how due to the translational invariance, the second order term vanishes in the 
homogeneous system. 

4.1 Translationally Invariant lattice 

For the translational invariant lattice we use the quasiparticle transformation 
given by Eqs. (26) and (27). Under this transformation the self consistent 
equations can be written, generalizing the previous version, as: 



n = (V| 2 + y]TK| 2 ^-2j, ( 5 °) 

ki 2 = m+i^i+^i, (5i) 
K|2 = m±miz^i, (52) 

N = 7|z| 2 + ^|^| 2 . (53) 

Q 

Here the phonon excitation spectrum is given by: 

huj q - y/K(q)[2\z\*V + K(q)], (54) 

and K(q) is given by Eq. (35). Again we omit the subscript in the amplitudes 
Zi because they have the same value at all lattice sites. Notice also that due to 
translational invariance, the c q coefficients vanish. 

In the homogeneous system, the form of the HFB-Popov equations for the 
quasiparticle amplitudes and energies is very close to the standard Bogoliubov 
form. We do, however, have to replace n = N/I by the condensate amplitude 
| z | 2 which must take into account the depletion of the condensate. We solve 
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for the condensate amplitude as a function of the external parameters J, V, N 
and I by inserting Eq. (52) in Eq. (53). Once \z\ 2 is determined, we use it to 
calculate the other expressions. 

In Fig. 1 we compare the number fluctuations on a lattice site, the con- 
densate fraction and the total and second order superfluid fraction determined 
from the exact solution of the Bose-Hubbard Hamiltonian to the self consistent 
HFB-Popov predictions as a function of the ratio V e ff = V/ J. The systems 
used for the comparisons have three wells, 7 = 3, and commensurate filling 
factors n = 5, 10, 20 and 50. Similar results for the incommensurate case with 
N = 16,31,61, 151 are shown in Fig. 2. We were restricted to consider only 
three wells due to computational limitations. The size of the matrix needed in 
the exact solution for N atoms and I wells scales as ^T7^W • However, if the 
HFB-Popov approach works well for these small systems we expect it to provide 
a good description of the larger systems prepared in the lab. 

Because the second order term of the superfluid fraction (second term of 
Eq.14) vanishes in the HFB-Popov approach (see Eq. 48), we only expect the 
self consistent HFB-Popov theory to give a good description of the superfluid 
fraction in the region where the second order term is extremely small, provided 
it predicts accurately the first order term. This is exactly what is observed in the 
plots. When the second order term starts to grow, typically above 0.5V°'^j, the 
HFB-Popov theory starts to fail. An estimate of Vffi is shown by a vertical line 
in some of the figures. This was obtained by using the second order perturbative 
approach presented in Ref . [4] . With increasing filling factor the critical value 
is shifted towards larger values of the interaction strength, and the region in 
which the HFB-Popov theory is accurate gets larger. It is interesting to note 
that the number fluctuations predicted by the theory are accurate in a greater 
range than the other physical quantities shown. Its predictions of squeezing 
agree very well with the exact solutions right up to the point where the number 
fluctuations become less than unity. 

For the cases with non-commensurate fillings depicted in Fig. 2, the agree- 
ment is significantly better for all quantities. This is not surprising because 
when the filling is not commensurate there is always a superfluid present and 
the Mott transition doesn't occur. As can be seen in the plots for these cases 
the second order term is always very small. 

4.2 Inhomogeneous lattice 

In this section we consider the experimentally relevant case when there is an 
external magnetic confinement in addition to the lattice potential. In this situ- 
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Figure 1 : Comparisons of the exact solution (dotted line) and HFB-Popov (solid 
line) as a function of V e // = V/J, for a system with 1 = 3 and filling factors 
n D = 5, 10, 20, and 50. Top: number fluctuations, middle: condensate fraction, 
bottom: superfluid fraction / s . The exact second order term (dashed line) of 

the superfluid fraction, / s is also shown in these plots. The vertical line shown 
in some plots is an estimation of V^f ■ 



ation, the self consistent HFB-Popov equations take the form: 



flUJ q U q i +C q Zi 



-TVjJqVl 



hi 
N 



= (2V(\z l \ 2 +n i ) - M + Sli 2 )ul- J(u q i+1 +uU) - Vz 2 lV %55) 
= (2V(\zf+h t ) - M + nt 2 )v q ~ J(v q i+1 +vU) - Vz?u$6) 



-J(z i+1 +z i -. 1 )+ (V(\z t \ 2 



2h, 



i 



* q 



ZiV q ) 



(57) 
(58) 

(59) 

(60) 



where VI = ^muj 2 a 2 , with m the mass of the atoms, u> the trap frequency, and a 
the lattice spacing. The site indices i are chosen such that i — corresponds to 
the center of the trap. Again the c q are introduced to ensure the orthogonality of 
the excitations to the condensate [23] . We have solved the HFB-Popov equations 
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Figure 2: Comparisons of the exact solution (dotted line) and HFB-Popov 
(solid line) for a system with 1 = 3 and non commensurate filling factors 
N = 16,31,61,151 as a function of V e ff = V/J. Top: number fluctuations, 
middle: condensate fraction, bottom: superfluid fraction / s . In these plots the 
exact second order term of the superfluid fraction is also shown with a dashed 
line. 



for this system by an iterative procedure, similar to the one followed in Rcf. 
[24]. Each cycle of the iteration consists of two steps. In the first step we 
solve Eq. (57) subject to the constraint Eq. (59) by using the hi obtained in 
the previous cycle. This generates new values for the In the second step 
we solve for in Eqs. (55) using the hi from the previous cycle and the 

newly generated Zi. The {u?, v\) are used then to update hi. Because the HFB- 
Popov is gapless, it is possible to keep the orthogonality of the excitations to 
the condensate by solving Eqs. (55) with the c q set to zero but removing in each 
cycle the projection of the calculated {uf,vf} amplitudes onto the condensate. 
Convergence is reached when the change in J2i l™i| 2 horn one cycle to the next 
is smaller than a specified tolerance. 

The parameters chosen for the numerical calculations were il = 0.0015-Er, 
with Er the one photon recoil energy, which for the case of a Rubidium con- 
densate corresponds to a trap frequency of approximately 90 Hz. We used a 
total number of 1000 atoms, N = 1000, and set VN = I.OEr. J was varied to 
achieve a range of V e ff = V/J between 0.01 and 312. The range was chosen 
based on a local mean field approach [4], which for our parameters estimates 
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the transition region between V e ff ~ 640 (at the center where the local filling 
factor is approximately 80) and V e ff ~ 12 (at the wings). 

The results of the numerical calculations are summarized in Figs. 3 to 8. In 
Fig. 3 we plot the evolution of the density profile (black boxes), the condensate 
population (triangles) and the on-site depletion (empty diamonds) as V e ff is 
increased. In the plots we also show, for comparison purposes, the ground state 
density profile for J = (empty boxes). This has the advantage that can be 
calculated exactly from the Hamiltonian. In general we observe the reduction 
of the condensate population and thus the increment of the depletion as the 
interaction strength is increased. When the system is in the superfluid regime 
most of the atoms are in the condensate but as J is decreased the depletion of 
the condensate becomes very important. 

For the chosen parameters, the density profile has a parabolic shape reflecting 
the confining potential. By comparing the evolution of the density as J is 
decreased with the exact solution at J — 0, we can crudely estimate the validity 
of the HFB-Popov calculations. The density evolves from a Gaussian type (see 
plots for V e ff — 0.01 and 0.09) with smooth edges towards a Thomas- Fermi 
profile with sharp edges adjusting its shape to the J = profile. We can 
appreciate that around V e ff — 3 both profiles are almost equal. For lower 
values of J the HFB-Popov density starts to differ from the J = one, even 
though the system is closer to the J = limit. We can say that beyond this 
point higher order correlations, neglected by the theory, begin to be important. 
The departure of the HFB-Popov density profile from the J = one as J 
is decreased begins at the edges (see panel corresponding to Veff = 11 and 
100). This is something expected if we look at the on-site depiction. For such 
values of V e ff the local depletion at the wings corresponds to a considerable 
percentage of the condensate populations, and thus the validity of the HFB- 
Popov assumptions starts to be dubious. The homogeneous results shown in the 
previous section corroborate our present statements for the confined system. For 
the smallest filling factor (see Fig.l) the differences between the homogeneous 
HFB-Popov calculations and the exact solutions become important for values 
of Veff greater than 20. For higher values of V e ff, see plot for V e ff — 312, 
the HFB-Popov density predictions differs from the J = solution even at the 
central wells. At this point the failure of the method is clear and a fully quantal 
method is required. 

The HFB-Popov quasiparticle spectrum is shown in Fig. 4. It can be observed 
how the lower energy eigenvalues evolve from a linear non degenerated spectrum 
to an almost degenerated one as J is decreased. It is worth it to mention that 
the small energy difference between the ground and first excited states for high 
values of V e ff makes the numerical solution very unstable in the sense that it 
was very easy to jump to an excited state when solving for the condensate wave 
function. The decrement in the energy spacing predicted by the HFB-Popov 
theory as the system approaches the transition is very useful to keep in mind 
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for the experimental realization of the Mott transition. As the optical lattice 
depth is ramped up the adiabaticity criteria is harder to fulfill. 

In Fig. 5 we plot the results for the number fluctuations found numeri- 
cally using the inhomogeneous HFB-Popov approach. The number fluctuations 
profile reflects the condensate profile. We also show the number fluctuations 
evaluated by using a local density approximation (empty boxes). The latter 
was calculated by substituting in the number fluctuations expression (Eq. 49) 
the {u q ,v q } amplitudes found for the homogeneous system (Eqs. 51 and 52), 
but replacing the condensate density in each lattice site by the one found nu- 
merically for the trapped system (see Fig. 3). The complete agreement between 
the two approaches justifies the validity of the local density approximation for 
the estimations of local quantities in confined systems. Based on this agreement 
and the results for the homogeneous system shown in the previous section, we 
expect that the inhomogeneous HFB-Popov results for squeezing also agrees 
with the exact solution right up to the transition. 



4.3 Quasi-momentum distribution in the inhomogeneous 



The quasi-momentum distribution of the atoms released from the lattice is im- 
portant because it is one of the most easily accessible quantities to the experi- 
ments. The quasi-momentum distribution function n q is defined as [10] 



where the quasi-momentum q can assume discrete values which are integer mul- 
tiples of a is the lattice spacing. In Fig. 6 we present the quasi-momentum 
distribution for the same parameters used in the previous section. The distribu- 
tion for the two lowest values of V e f / corresponds to the one that characterizes 
an uncorrelated superfluid phase with a narrow peak at small quasi-momenta. 
As the hopping rate is decreased we observe that the sharpness of the central 
peak decreases and the distribution extends towards large quasi-momenta. It is 
interesting to note the appearance of a small peak between q = 0.5 and 1 which 
is most noticeable for the V e fj — 3 case. This agrees with the results found 
in [25] where they solve numerically the Bose-Hubbard Hamiltonian by using 
Monte Carlo simulations. We attribute the origin of the small peak to the deple- 
tion of the condensate at the wings. For the parameters when the small peak is 
present, the most important contribution to the quasi-momentum distribution 
still comes from the condensate atoms. The step function like shape of the con- 
densate profile causes an oscillatory |sin(a;)/a;| shape of the quasi-momentum 
distribution. As the lattice depth is increased the hopping becomes energeti- 
cally costly, the long-range order starts to decrease and the Fourier spectrum 
becomes broader. 



lattice 




(61) 
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In Fig. 7 we plot the first order on site supcrfluid fraction f si which was 
defined in Eq.(48). The curves corresponding to V e ff = 0.01 — 11, which are 
in the regime where the HFB-Popov is expected to be valid, depict how as 
V e f f is increased the superfluid profile decreases faster at the wings and at the 
center but no major change is observed in the middle section. The evolution 
of the on-site superfluidity as the interaction strength is increased, exhibiting 
a domain localized decrement instead of a global one, is in agreement with the 
development of uncompressible regions surrounded by superfluid rings predicted 
for trapped systems [26] as the transition is approached. 

Fig. 8 shows the first order and total supcrfluid fraction and also the second 
order superfluid fraction as a function of V e ff. Different from the translation- 
ally invariant case, the second order contribution calculated in the HFB-Popov 
approach doesn't vanish for the inhomogcncous system. The rapid decrement of 
the superfluid fraction observed after Log(V e ff) <~ 1.2 is a signature that the 
system is entering a highly correlated regime. Beyond this point higher order 
correlations neglected in the HFB-Popov approach become crucial and a more 
accurate approach is required. 

The Mott transition is a quantum phase transition and as all critical phe- 
nomena its behavior depends strongly on the dimensionality of the system. In 
the present analysis, due to computational limitations, we considered one di- 
mensional systems. Experimentally, the Mott transition has been achieved [2] 
in a 3 dimensional lattice with filling factors between 1 and 3. Even though the 
HFB-Popov approach fails to describe the strong coupling regime for the one 
dimensional systems we considered in the present paper, we showed how the 
method is incredibly powerful in describing most of its characteristic features as 
they are driven from the supcrfluid regime towards the transition. We expect 
the HFB-Popov method to give a better description of the transition as the 
dimensionality of the system is increased and therefore to be a good model in 
an experimental situation. 

As shown in previous studies [27], [28] the Mott transition in a d-dimensional 
homogeneous system has two different critical behaviors: one (d+1) XY- like, 
for systems with fixed integer density as the interaction strength is changed, 
and one mean field-like exhibited when the transition is induced by changing 
the density. Different from the homogeneous case where the Mott transition 
is characterized by the global offset of the superfluidity, for confined systems, 
commensuration is only well defined locally. The inhomogeneity introduced by 
the confined potential allows the existence of extended Mott domains (above a 
critical interaction strength) surrounded by superfluid ones [26], thus the total 
superfluid fraction doesn't vanishes in the Mott regime. This issue, together 
with the fact that the finite length scale introduced by the trap suppresses the 
long wave fluctuations which are responsible for destroying the mean field [23] 1 , 
make us believe the critical behavior in confined systems to be more mean-field 
like. Because the critical dimension for the latter type of transition is two [27], 

1 One obvious consequence of this is that BEC is possible in one and two dimensions in a 
trap whereas in the homogeneous, thermodynamic limit it can not occur in fewer than three 
dimensions 
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[28], we expect that for trapped systems in d = 3, the range of validity of the 
HFB-Popov extends closer to the transition. 

5 Summary 

We have developed in this article a Bogoliubov method for describing the ap- 
proach of a condensate loaded in an optical lattice towards the Mott transition. 
We have shown that this method can be used to predict the relevant physical 
quantities over a useful range. We have also shown how it gives a powerful 
insight into the way quantum depletion reduces the long range order and the 
supcrfluid fraction. 
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Figure 3: Condensate density (triangles), total density (filled boxes) and local 
depletion (empty diamonds) as a function of the lattice site for different values of 
V e ff. Although these quantities are defined only at the discrete lattice sites we 
join them to help visualization. The empty boxes represent the exact solution 
for the case J=0. 
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Figure 4: Quasiparticle spectrum in ascending order predicted by the HFB- 
Popov theory for different values of V e ff. Empty diamonds (V e /f — 0.01), stars 
(V e ff = 0.09), crosses (V e ff = 3), filled diamonds (V e ff = 11), empty boxes 
(V ef f = 100) and polygons (V ef f = 312) 
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Figure 5: Number fluctuations in the self consistent HFB-Popov approach as a 
function of lattice site for V e ff = 0.01 (boxes), V e ff = 0.09 (crosses), V e // = 3 
(circles), 1 ^// = H (triangles), V e // = 100 (stars) and V e ff = 312 (diamonds). 
The maximum value reached by the profile decreases as V e ff is increased. The 
empty boxes shown for each of the curves correspond to the number fluctua- 
tions predicted by the homogeneous HFB-Popov model using a local density 
approximation. 
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Figure 7: First order on-site superfluid fraction as a function of the lattice site for 
different values of V e ff. Filled boxes: V e ff. =0.01, empty boxes: 14//. = 0.09, 
empty diamonds: 14//. = 3, stars: 14//. = 11, crosses: 14//. = 100 and 
triangles: 14//. = 312. 
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Figure 8: Top panel: First order (boxes) and total (stars) superfluid fraction 
as a function of V e ff. ■ Bottom panel : second order superfluid fraction as a 
function of V e f /. All these quantities are calculated in the self consistent HFB- 
Popov approach. 
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